setwd(out_path)
if(!dir.exists("Outcomes")){
	stop("Run model estimation first")
}else{
	setwd("Outcomes")
	if(!dir.exists("post")){
		dir.create("post")
	}else{
		warning("Overwriting previous run, if any")
	}
}
setwd("post")

################
#Predicting policy passage from partisan support
################
minmax_seq<-function(numvec){
	.out<-seq(from=min(numvec),
		to=max(numvec),
		by=(diff(range(numvec)/10))
	)
	return(.out)
}

############
#Calculate expected value, with confidence intervals, for GLM
############
glm_expval<-function(glm.model,sim.data){
	#Calculate expected value, with confidence intervals, for GLM
	
	#Coef matrix
	coefMCraw<-coef(glm.model)
	# coefMC2rm<-grep("factor(GovtControl)",names(coefMCraw),value=T,fixed=T)
	# coefMC2rm<-coefMC2rm[!grepl("Dem",coefMC2rm,fixed=T)]
	
	# coefMC<-coefMCraw[setdiff(names(coefMCraw),coefMC2rm)]
	coefMC<-coefMCraw
	coefMC.mat<-matrix(coefMC,nrow=length(coefMC),dimnames=list(names(coefMC)))
	
	#X matrix
	X.mat<-as.matrix(sim.data[,names(coefMC)[2:(length(coefMC))]])
	X.mat<-cbind("(Intercept)"=1,X.mat)
	class(X.mat)<-"numeric"
	
	#Calculate linear predictor
	stopifnot(dim(X.mat)[2]==dim(coefMC.mat)[1])
	stopifnot(all(colnames(X.mat)==rownames(coefMC.mat)))
	Xb<-X.mat%*%coefMC.mat
	
	# #Calculate cutpoints
	# cutpointsMC.full<-rep.full.op$zeta
	# cutpointsMC<-matrix(cutpointsMC.full,nrow=1)
	
	#Get confidence intervals
	SEs<-as.matrix(sqrt(diag(X.mat%*%vcov(glm.model)[names(coefMC),names(coefMC)]%*%t(X.mat)))*1.96)
	
	#Calculate fits
	lfit<-Xb
	lfit.hi<-lfit+SEs
	lfit.lo<-lfit-SEs
	fit<-exp(lfit)/(1+exp(lfit))
	fit.hi<-exp(lfit.hi)/(1+exp(lfit.hi))
	fit.lo<-exp(lfit.lo)/(1+exp(lfit.lo))
	pred.df<-data.frame(fit=fit,fit.hi=fit.hi,fit.lo=fit.lo)
	return(pred.df)
	
}

#Fake Data for Graphs (Democratic Support)
sim.data<-data.frame(b.scaled.intgrp =  mean(aic_clean$b.scaled.intgrp) ,
					 m.scaled.intgrp =  mean(aic_clean$m.scaled.intgrp) ,
                     pred50_sw = mean(aic_clean$pred50_sw) , 
                     pred90_sw = mean(aic_clean$pred90_sw) ,
                     rescaled_ideol_code = mean(aic_clean$rescaled_ideol_code, na.rm = TRUE) ,
                     fp.policies = 1,
					 econ.policies = 0,
					 rep = 0,
					 dem = minmax_seq(aic_clean$dem))
ohat_90_ep <- glm_expval(full.mv.glm,sim.data)
ohat_aff <- cbind(as.data.frame(ohat_90_ep), party=sim.data$dem, id="Democratic Party Support")


#Fake Data for Graphs (Republican Support)
sim.data<-data.frame(b.scaled.intgrp = mean(aic_clean$b.scaled.intgrp) ,
					 m.scaled.intgrp = mean(aic_clean$m.scaled.intgrp),
                     pred50_sw = mean(aic_clean$pred50_sw) , 
                     pred90_sw = mean(aic_clean$pred90_sw) ,
                     rescaled_ideol_code = mean(aic_clean$rescaled_ideol_code, na.rm = TRUE) ,
                     fp.policies = 1,
					 econ.policies = 0,
					 rep = minmax_seq(aic_clean$rep),
					 dem = 0)
ohat_50_ep <- glm_expval(full.mv.glm,sim.data)
ohat_mc <- cbind(as.data.frame(ohat_50_ep), party=sim.data$rep, id="Republican Party Support")

#Now Let's melt the dfs
ohat.full <- rbind(ohat_mc, ohat_aff)
op.plot <- ggplot(ohat.full)

#FP Graph
post.passage.partisan <- op.plot +
  geom_point(aes(x=party,y=fit,group=id,shape=id)) +
  geom_line(aes(x=party,y=fit,group=id,linetype=id)) + 
  geom_ribbon(aes(x=party,ymin=fit.lo,ymax=fit.hi,fill=id),alpha=.2)+
  ggtitle('Predicted probability of policy adoption',sub='Based on partisan support') +
  xlab(" Support for Proposed Policy") +
  ylab(" Probability of Policy Adoption") +
  theme_bw()

ggsave("post_passage_partisan.tiff",post.passage.partisan)

if(F){
	#Validate that it predictions are non-linear
	pdf("post_passage_partisan_VAL50.pdf")
	plot(diff(ohat_50_ep$fit))
	title("Difference plot of predicted values, middle class")
	dev.off()
	
	#Validate that it is non-linear
	pdf("post_passage_partisan_VAL90.pdf")
	plot(diff(ohat_90_ep$fit))
	title("Difference plot of predicted values, affluent")
	dev.off()
}

######################
#Predicting policy adoption from interest group support
######################

#Fake Data for Graphs (Business Opinion)
sim.data<-data.frame(b.scaled.intgrp =  minmax_seq(aic_clean$b.scaled.intgrp) ,
					 m.scaled.intgrp =  mean(aic_clean$m.scaled.intgrp) ,
                     pred50_sw = mean(aic_clean$pred50_sw) , 
                     pred90_sw = mean(aic_clean$pred90_sw) ,
                     rescaled_ideol_code = mean(aic_clean$rescaled_ideol_code, na.rm = TRUE) ,
                     fp.policies = 1,
					 econ.policies = 0,
					 rep = median(aic_clean$rep),
					 dem = median(aic_clean$dem))
ohat_90_ep <- glm_expval(full.mv.glm,sim.data)
ohat_aff <- cbind(as.data.frame(ohat_90_ep), scaled_ig=sim.data$b.scaled.intgrp, id="Business Interest Group Alignment")


#Fake Data for Graphs (Middle Class Opinion)
sim.data<-data.frame(b.scaled.intgrp = mean(aic_clean$b.scaled.intgrp) ,
					 m.scaled.intgrp = minmax_seq(aic_clean$m.scaled.intgrp),
                     pred50_sw = mean(aic_clean$pred50_sw) , 
                     pred90_sw = mean(aic_clean$pred90_sw) ,
                     rescaled_ideol_code = mean(aic_clean$rescaled_ideol_code, na.rm = TRUE) ,
                     fp.policies = 1,
					 econ.policies = 0,
					 rep = median(aic_clean$rep),
					 dem = median(aic_clean$dem))
ohat_50_ep <- glm_expval(full.mv.glm,sim.data)
ohat_mc <- cbind(as.data.frame(ohat_50_ep), scaled_ig=sim.data$m.scaled.intgrp, id="Advocacy Interest Group Alignment")

#Now Let's melt the dfs
ohat.full <- rbind(ohat_mc, ohat_aff)
op.plot <- ggplot(ohat.full)

#FP Graph
post.passage.interestgroups <- op.plot +
  geom_point(aes(x=scaled_ig,y=fit,group=id,shape=id)) +
  geom_line(aes(x=scaled_ig,y=fit,group=id,linetype=id)) + 
  geom_ribbon(aes(x=scaled_ig,ymin=fit.lo,ymax=fit.hi,fill=id),alpha=.2) + 
  ggtitle('Predicted probability of policy adoption',sub='Based on interest group support') +
  xlab(" Support for Proposed Policy") +
  ylab(" Probability of Policy Adoption") +
  theme_bw()

ggsave("post_passage_interestgroups.tiff", post.passage.interestgroups)

if(F){
	#Validate that it is non-linear
	pdf("post_passage_interestgroups_VAL50.pdf")
	plot(diff(ohat_50_ep$fit))
	title("Difference plot of predicted values, middle class")
	dev.off()
	
	#Validate that it is non-linear
	pdf("post_passage_interestgroups_VAL90.pdf")
	plot(diff(ohat_90_ep$fit))
	title("Difference plot of predicted values, affluent")
	dev.off()
}

################
#Get odds ratios from full models
################
mod.data<-function(mod){
	mod.df<-data.frame(summary(mod)$coefficients)[-1,1:2] #drops intercept & factors
	mod.df$varRaw<-rownames(mod.df)
	return(mod.df)
}

relabel_vars<-function(invec){
	varnames<-c(MidClassPref="pred50_sw",AffluentPref="pred90_sw",
				BusinessPref="b.scaled.intgrp",AdvGroupPref="m.scaled.intgrp",
				IdeolCode="rescaled_ideol_code",
				EconPolicy="econ.policies",ForPolicy="fp.policies",
				RepSupport="rep",DemSupport="dem"
				)
	stopifnot(all(invec %in% varnames))
	newvarloc<-match(invec,varnames)
	newvars<-names(varnames)[newvarloc]
	stopifnot(length(newvars)==length(invec))
	newvars<-factor(newvars,
				levels=rev(names(varnames)[names(varnames) %in% newvars]),
				labels=rev(names(varnames)[names(varnames) %in% newvars])
				)
	return(newvars)
}
passage_data<-mod.data(full.mv.glm)
passage_data$var<-relabel_vars(passage_data$varRaw)
passage_data$EstLo<-(passage_data$Estimate-(1.96*passage_data$Std..Error))
passage_data$EstHi<-(passage_data$Estimate+(1.96*passage_data$Std..Error))

# passage_data$Estimate<-exp(passage_data$Estimate)
# passage_data$EstLo<-exp(passage_data$Estimate-(1.96*passage_data$Std..Error))
# passage_data$EstHi<-exp(passage_data$Estimate+(1.96*passage_data$Std..Error))


#Plot coefficients
passage.coef<-ggplot(passage_data)+
	geom_pointrange(aes(y=Estimate,ymin=EstLo,ymax=EstHi,
			x=var),size=0.2)+
	geom_hline(aes(yintercept=0),size=0.2)+
	theme_bw()+ylab("Logit coefficient")+xlab("Feature")+
	ggtitle("Logit parameters for policy passage",subtitle="With 95% confidence intervals")+
	coord_flip()

ggsave("post_passage_coef.tiff",passage.coef)

################
#Get coeficients from models, by issue
################
mod.data<-function(mod){
	mod.df<-data.frame(summary(mod)$coefficients)[-1,1:2] #drops intercept
	mod.df$varRaw<-rownames(mod.df)
	return(mod.df)
}

passage_econ_data<-mod.data(ec.full.mv.glm)
passage_econ_data$var<-relabel_vars(passage_econ_data$varRaw)
passage_econ_data$IssueArea<-"Economic"

passage_fp_data<-mod.data(fp.full.mv.glm)
passage_fp_data$var<-relabel_vars(passage_fp_data$varRaw)
passage_fp_data$IssueArea<-"ForeignPolicy"

passage_ne_data<-mod.data(ne.full.mv.glm)
passage_ne_data$var<-relabel_vars(passage_ne_data$varRaw)
passage_ne_data$IssueArea<-"Social"

passage_byissue_data<-rbind(passage_econ_data,passage_fp_data,passage_ne_data)
jitterspace<-sd(passage_byissue_data$Std..Error)/3

#Plot coefficients
passage.byissue.coef<-ggplot(passage_byissue_data)+
	geom_pointrange(aes(y=Estimate,ymin=(Estimate-(Std..Error*1.96)),ymax=(Estimate+(Std..Error*1.96)),
			x=var,linetype=IssueArea,shape=IssueArea),position=position_dodge(jitterspace),size=0.2)+
	geom_hline(aes(yintercept=0),size=0.2)+
	theme_bw()+ylab("Parameter estimate")+xlab("Feature")+
	ggtitle("Paramter estimates for policy passage, by issue area",subtitle="With 95% confidence intervals (Maximum likelihood estimation)")+
	coord_flip()

ggsave("post_passage_byissue_coef.tiff",passage.byissue.coef)

setwd(script_path)